Albendazole repurposing on VEGFR-2 for possible anticancer application: In-silico analysis

Drug repurposing is the finding new activity of the existing drug. Recently, Albendazole’s well-known antihelmintic has got the attention of an anticancer drug. Plausible evidence of the interaction of Albendazole with one of the types of tyrosine kinase protein receptor, vascular endothelial growth factor receptor-2 (VEGFR-2) is still not well understood. Inhibition of the VEGFR-2 receptor can prevent tumor growth. The current study investigated the interaction of Albendazole with VEGFR-2.It was found that the said interaction exhibited potent binding energy ΔG = -7.12 kcal/mol, inhibitory concentration (Ki) = 6.04 μM, and as positive control comparison with standard drug (42Q1170A) showed ΔG = -12.35 kcal/mol and Ki = 881 μM. The key residue Asp1046 was formed involved hydrogen bonding with Albendazole. The molecular dynamics simulation study revealed the stable trajectory of the VEGFR-2 receptor with Albendazole bound complex having significant high free energy of binding as calculated from Molecular Mechanics Generalized Born and Surface Area study ΔG = -42.07±2.4 kcal/mol. The binding energy is significantly high for greater stability of the complex. Principal component analysis of molecular docking trajectories exhibited ordered motion at higher modes, implying a high degree of VEGFR-2 and Albendazole complex stability as seen with the standard drug 42Q. Therefore, the current work suggests the role of Albendazole as a potent angiogenesis inhibitor as ascertained by its potential interaction with VEGFR-2. The findings of research will aid in the future development of Albendazole in anticancer therapy.


Introduction
Cancer is most deadly disease in the world. Cancer treatment is the most difficult task for researchers. As a result, researchers devote time and resources to discovering and developing new therapies for various cancers. Nonetheless, cancer occurrence is anticipated to rise by even more over 50% in the next years, which is regrettable. Ovarian cancer is the seventh leading cause of mortality and morbidity in women globally. There are many types of ovarian cancer, both common and uncommon, such as epithelial ovarian cancer, germ cell tumour, primary peritoneal cancer, fallopian tube cancer, stromal cell tumour, and ovarian cancer [1]. In general, cancer is treated using traditional methods that target the DNA of all cells; however, one important disadvantage of this treatment is that it also affects healthy cells. To address this issue, targeted drug treatment has been established to release a drug to a specific site of receptor while minimizing cell toxicity [2]. To design targeted treatments, scientists must first discover the genetic alterations that allow a cancer cell to mutated. A protein found in cancer cells may become promising approach to target tumor instead of healthy cell [3]. In the treatment of cancer, a broad variety of targeted therapies have been tried. Various drug administration techniques, such as hormone treatments, signal transduction inhibitors, gene expression modulators, apoptosis inducers, angiogenesis inhibitors, immunotherapies, and toxin delivery vehicles, have been shown to have stronger therapeutic action [3]. Angiogenesis inhibitor is one of the successful strategies from this group. Angiogenesis is required for tissue formation, regeneration, and reproduction and is a basic process in the development of a wide range of clinical diseases, including cancer [4]. Angiogenesis, the process by which new blood vessels enter tumour masses, giving oxygen and nutrients to help tumour development and spread, is also important in the transformation of a benign tumour into a malignant one [5]. VEGF has been recognised as an important controller of blood vessel formation and a major facilitator of tumor angiogenesis since its finding in 1983. VEGF is basically falls under growth factor and platelets derived growth factor family that responsible for vasculogenesis and angiogenisis. Overexpression of VEGF in serum revealed an increase in tumour vascular density, which correlates with the degree of malignancy. According to the findings, VEGF levels in ovarian cancer is increased by more than 70% and also in various cancer types. VEGF ordinarily mediates its actions via binding to VEGF receptor present on endothelial cells and by direct acting on VEGFR on tumor cell. Endothelial cell growth and migration can also be aided by VEGF. Endothelial cells are abundant in cancerous tissue [4,6]. VEGF also controlled number of functions in the body like, invasion, endothelial cell progression, migration, cell permeability and dilation of vessels. So angiogenesis in cancer conditions can be controlled by blocking the tyrosin kinase signalling pathway of VEGF. When VEGF attaches to VEGFR, it causes it to alter conformation, dimerize, and phosphorylate tyrosine [7]. Generally, VEGF acting on both VEGFR1 and VEGFR2 receptor, but more intracellular signalling activation shown by VEGFR2 than VEGFR1. Consequently, VEGFR-2 receptor-associated signalling has emerged as a prospective therapeutic target for cancer therapy [8]. As a result, inhibiting or down-regulating VEGFR-2 signalling is a significant method for creating novel therapeutics for many human malignancies [9]. Drugs that have been licensed or evaluated for other illnesses but demonstrate surprising cytotoxicity against tumor sites, on the other hand, might be repurposed as chemotherapeutic options [10]. Because such drugs are already well-known due to their pharmacokinetics and dynamics features, the odds of failure are quite low, and the development stage becomes highly economical [11]. Furthermore, the majority of antihelmintic drugs that have already been approved are being repurposed as antitumor agents [12,13]. Notably, fenbendazole, a new tubulin-binding agent with anticancer properties, has recently been proposed as a potential anticancer agent [10]. Mebendazole could also be used for tumor cell inhibition due to its pharmacokinetic features and microtubule depolymerization action [14].
The effectiveness of albendazole as anti-helminthic drug and other benzimidazole class of drug has been related mostly due to its tubulin binding capability, which results in parasite depolymerization and cell death. Because of their capacity to disrupt microtubule formation, benzimidazoles have been studied for their anti-tumor efficacy in different types of cancer [15]. Albendazole (ABZ) is an approved anthelmintic drug has an anticancer activity and also albendazole shown to have a wide effect on paclitaxel-resistant carcinoma cells [16]. Albendazole demonstrated numerous mechanisms of action in order to show its antihelmintic and anti-cancer potential Albendazole inhibits tubulin polymerization, which is one of the key action pathways. The suppression of tubulin production is directly connected to angiogenesis and vasculogenesis. This binding affects the cell cycle, which leads to cell death by apoptosis. Angiogenesis (the formation of new blood cells) is prevented by inhibiting VEGF receptor overexpression [17].
In this study, mainly focused on efficacy of albendazole as an anticancer agent by studied different pathways. Also studied probability of activity by using PASS prediction, finding showed good activity against carcinogeneic cell. Molecular docking, MD simulations were run to validate stable and converged complexes and learn more about the drug's thermodynamic characteristics at the specific receptor. The binding energy of ligands to protein molecules determined by MMGBSA (Molecular Mechanics Generalized Born and Surface Area) calculations and it showed significant affinity towards 3VHE target protein. Free surface energy of ABZ bound VEGFR2 receptor analyzed by using PCA (Principal component analysis).

Predictions of molecular targets
By entering smiles of desired drugs, the SwissTargetPrediction service was used to investigate specific targets gives an idea about their physiological side effects and cross reactivity for anticancer efficacy. A three-dimensional & two-dimensional similarity index with known ligands is used to determine the probable targets of the provided chemical [20].

Prediction of biological activity spectra for substances
Prediction of Activity Spectra for Substances (PASS) is a free web service that predicts the biological activities of drug-like chemicals (http://www.way2drug.com/passonline). Based on multilevel neighbourhoods of atom (MNA) descriptors and a modified Bayesian algorithm containing more than 300,000 organic compounds with more than 6,825 distinct biological activities, it also predicts a structure-activity relationship (SAR). Thus, PASS may be applied to improve and well-target chemical production and biological testing. The anticancer activity of the ABZ and standard drug 42Q was anticipated, and the findings are shown in Table 1. A substance's estimated activity is forecasted as probable activity (Pa) and probable inactivity (Pi). Only compounds with a Pa greater than Pi are considered suitable for a certain therapeutic activity [21].

Molecular docking
This research aims to shed insight on the binding properties of 3VHE protein using geometrically optimized chemicals like ABZ and standard 42Q. Before being employed in molecular interaction studies, proteins were checked for missing side-chain residues using the open Molecular Mechanics (MM) simulation tool (https://openmm.org/). Autodock v 4.2.6 was used for molecular docking investigations. The binding cavity was established using the co- https://doi.org/10.1371/journal.pone.0287198.g001 crystallized X-ray structure of VEGFR-2 protein from the RCSB PDB. Within three spaces of the co-crystallized ligand, the residue locations were computed. The energy was minimized using the steepest descent and conjugate gradient techniques after cavity selection. Finally, the receptor and target molecules were stored in pdbqt format after combining the nonpolar hydrogens. Grid boxes with a spacing of 0.3 were made. To get the minimum free energy of binding (G), docking studies of the protein-ligand complex were performed utilizing the Lamarckian Genetic Algorithm (LGA). Three replicates of molecular modelling experiments were carried out, with default parameters like number of solutions (50 in each), number of evaluations (2500000), maximum number of generations (2700) etc. After docking, the Root Mean Square Deviation clustering maps were created by re-clustering with clustering tolerances of 0.25, 0.5, and 1 in order to obtain the best cluster with the lowest energy score and a high number of populations, respectively.

Molecular dynamics simulation
The Desmond software Vs 2020.1 from Schrödinger was used to do MD simulations on docked complexes for ABZ and 42Q with protein (PDB ID: 3VHE). Triplicate sampling was performed with the same condition for each MD run to achieve better results. This system employed the OPLS-2005 force field [22][23][24] and an explicit solvent model with simple point charger (SPC) water molecules [25]. To neutralise the charge, Na+ ions were added. NaCl solution (0.15M) was used to simulate the physiological environment. To retrain over the proteinligand complex, the system was first equilibrated using NVT ensemble for 200 ps. This system was then exposed to a 12 ps run for equilibration and minimization using the NPT ensemble. The NPT ensemble was built using the Nose-Hoover chain coupling method [26]. All simulations were performed at a temperature of 27˚C, with a relaxation time of 1.0 ps and a pressure of 1 bar. A time step of 2 fs was chosen. For pressure regulation, the Martyna-Tuckerman-Klein chain coupling scheme [27] barostat technique with a relaxation duration of 2 ps was utilised. The particle mesh Ewald method [28] was used to determine long-range electrostatic interactions, with the Coulomb interaction radius set at 9 cm. For each trajectory, the bonded forces were estimated using the RESPA integrator with a time step of 2 fs. To measure the Where, ΔG bind = binding free energy, G complex = free energy of the complex, G protein = free energy of the target protein, and G ligand = free energy of the ligand. The MMGBSA outcome trajectories were analyzed further for post dynamics structure modifications.

Principal component analysis (PCA)
During a 200 ns simulation of the human VEGFR2 kinase domain (PDB ID: 3VHE) complexed with ABZ, PCA analysis was used to recover the global movements of the trajectories. A covariance matrix was generated as described to calculate the PCA. Conformational analysis of Albendazole in-bound complex was done by using 10 appropriate modes of conformation. The main component was calculated as trajectories, and a comparison of the first highest mode (PC1) versus the second highest mode (PC2), followed by PC2 and PC3, and finally, PC9 and PC10 modes were investigated. Geo measures v0.8 were used to calculate the principal components based on eigenvectors from the correlation matrix on an ABZ-bound complex. The MD trajectory versus principal components of trajectory motion was recorded in a 2D plot using the Matplotlib python package using Geo measures, including a comprehensive library of g_sham [29].

Molecular target predictions
Molecular target studies were projected based on their similarities to known compounds to estimate their targets to conduct the molecular mechanisms behind a certain phenotypic or bioactivity and rationalize any side effects. The top 50 results of the closely related receptors in 2D/3D were shown as a pie chart based on Target, Popular Name, Uniprot ID, ChEMBL-ID, Target Class, Probability, and Reported actives (Fig 2). As can be seen, ABZ predicts 26% Kinase receptor, 18% family α G protein-coupled receptor, and 18% Enzymes while 42Q has 34% Kinases and 14% protease, and 6% Enzymes.
Overall, the findings demonstrated the efficacy of ABZ as an anticancer treatment that targets the MAP kinase p38 alpha, protein-tyrosine kinase erbB-2, nuclear factor-kappa kinase beta subunit, Serine/threonine-protein kinase, Tyrosine-protein kinase JAK3, Tyrosine-protein kinase JAK2, VEGFR-1.

Prediction of activity spectra for substances (PASS) prediction
PASS employed both ABZ and 42Q to forecast the biological spectrum. Out of a maximum likelihood score of 1, it calculates the probability of activity and inactivity against cancer and non-cancerous cells. With active coefficients of 0.348, 0.320, 0.248, and 0.212, ABZ showed considerable anticarcinogenic action against melanoma, solid tumors, lymphocytic leukemia, and endocrine cancer in the PASS filter (Table 1). 42Q was shown to be effective against melanoma (0.302) and solid tumors (0.348). Nonetheless, both ABZ and 42Q displayed significant activity on endocrine cancer; Platelet-derived growth factor receptor kinase, MAP3K5, Transcription factor STAT, Angiogenesis, etc. These findings point to a high likelihood of anticarcinogenic action. This bioactivity investigation confirms albendazole's importance in protecting human health against tumor development and inflammation.

Molecular docking
All the binding energy scores are calculated from the best cluster (95%) that falls within the lowest RMSD 0.25 Å. With the lowest binding energy, (G-7.12 kcal/mol) and inhibitory concentration, Ki (6.04 M), ABZ showed a considerable binding affinity for the human VEGFR2 kinase domain (PDB ID: 3VHE). During the interaction of the ligand ABZ, Asp1046 residues at the binding cavity of the protein were involved in conventional hydrogen bonding. Apart from hydrogen bonding, van der Waal's interactions by amino acid residues were also involved in non-bonded interaction with the ligand (Fig 3A and 3B left panel). While the standard ligand 42Q displayed a high binding affinity for 3VHE protein where free energy binds were recorded as -12.35 kcal/mol and inhibitory concentration (Ki) 881 pM, respectively. Glu885, Cys919, Asp1044 residues in the binding cavity of3VHE formed conventional hydrogen bonds (Fig 3B, right panel and 3C).

Molecular dynamics and simulation
The stability and convergence of the 3VHE+ABZ and 3VHE+42Q complexes were investigated using molecular dynamics and simulation (MD). When comparing the root mean square deviation (RMSD) data, a simulation of 200 ns revealed stable conformation. The RMSD of Cα-backbone of 3VHE bound to ABZ exhibited a deviation of 0.2 Å (Fig 4A (i)), while 42Q displayed a deviation of 0.64 Å. However, 3VHE protein without ligand (3VHE-Apo) exhibited a deviation from the beginning to the end of simulation 0.63 Å. Similar RMSD of ligand bound protein and apo-protein suggested least perturbations and stable conformations. Therefore, it can be suggested that ABZ and 42Q bound 3VHE is quite stable in complex due to the higher affinity of the ligand. Comparison of RMSD clustering having highest population of Apo protein and ABZ bound to 3VHE proteins displayed little conformational changes ( Fig  4A (ii)). At the binding cavity of ligand ABZ, couple of residue forming conventional hydrogen bonds viz. Glu917 and Cys919 where we have observed Glu917 pushed away (distance 2.1 Å) from the original position in ligand bound state while comparing with the position in apo protein where the orientation of the residue is much straighter (Fig 4A (ii), left panel, up). RMSD of the deviated residue of ligand bound state and apo state was measured to be RMSD 0.23 Å. The other residue Cys919 is pulled toward ligand ABZ(distance 1.1 Å) while in the apo protein Cys919 seemed to be straighter and unaltered (RMSD 0.21 Å) (Fig 4A (ii), left panel, up). On the other hand, at the C-terminal end residues Asp994-Lys997 of the apo protein conformed

PLOS ONE
randomly into a loop structure whereas, in ABZ bound 3VHE protein a partial helix formed at that position (RMSD 1.5Å) (Fig 4A (ii), left panel, down). Formation of helix structure made the C-terminal of 3VHE more stable than the loop orientation. The top three clusters having highest populations are also compared to see the positional variances of the ligand. Comparison of conformation of top two cluster 1 and cluster 2 displayed aligned ligands having small positional deviation (RMSD 0.25Å) (Fig 4A (ii), right panel, green-cyan, up). Whereas, comparison of cluster 1 and cluster 3 showed significant deviation of ligand twisted around the axis (Fig 4A (ii), right panel, green-magenta, down). Therefore, conformational stability is achieved in the untwisted form of the ligand at the binding cavity than the twisted conformation. Apo and 42Q bound top most clusters having highest population showed significant movements of residues to achieve high stability. In 42Q bound 3VHE (Fig 4A (ii), cyan) two key residues Cys919 and Asp1046 moved apart to form spatial arrangement with the ligand for hydrogen bond formation. While Glu885 pulled toward the ligand to attain a stable hydrogen bond formation where in apo protein the residue lie distantly apart (Fig 4A (ii)). Comparative analysis of top two cluster 1 and cluster 2 having highest population of 42Q bound 3VHE trajectories displayed least ligand mobility from its axis exhibiting very stable complexes (Fig 4A  (ii)). The plots for root mean square fluctuations (RMSF) displayed small spikes of fluctuation except at amino acid residue 60 (4.5 Å) and 130-140 (8 Å) in 3VHE-ABZ complex, while the rest of the residues less fluctuated during the entire 200 ns simulation (Fig 4B, dark pink) indicating the stable amino acid conformations during the simulation time. The high fluctuations in ABZ bound protein residues due to loop conformation. However, in 42Q bound protein showed significant fluctuation at 130-140 residues (4.3 Å) (Fig 1B, dark green). In contrast, RMSF values of the 3VHE-apo displayed significant fluctuations, and at 130-140 residues, the fluctuation was found to be 9 Å, signifying very flexible unstable residues ordered into a loop conformation (Fig 4B, orange). Therefore, from RMSF plots, it can be suggested that the structures of 3VHE were stable during simulation in ABZ and 42Q bound conformations. The number of hydrogen bonds between protein and ligand suggests the significant interaction and stability of the complex. A significant number of hydrogen bonds showed between 3VHE & ABZ and 42Q throughout the simulation time 200 ns (Fig 4C). Consistent numbers of hydrogen bonds are observed between 3VHE and ABZ (Average 2 numbers) and 42Q (Average 4 numbers) that might facilitate to conform into stable complex (Fig 4C). The gyration radius is a measurement of the protein's compactness. Here in this study, 3VHE Cα-backbone bound to ABZ displayed lowering of radius of gyration (Rg) from 20.3 Å to 20.0 Å (Fig 4D,  dark pink), and 3VHE Cα-backbone bound to 4Q2 displayed stable as well as lowering of peaks from 20.2 to 20.0 Å (Fig 4D). Similarly, 42Q bound protein displayed the lowering of Rg from 20.3 to 20.2 Å (Fig 4D). Lowering of RMSD indicates the compactness of protein ligand complex. On the other hand, apoprotein (3VHE-apo) displayed significant uneven fluctuations and was devoid of stability, indicating the less compact orientation of the protein ( Fig  4D, orange). After Rg analysis, similar patterns were identified in solvent accessible surface area (SASA) in both ligand-bound and unbound states. The unbound state of ligand 3VHE revealed a large surface area accessible to solvent (Fig 4E and 4F, black), but binding with 42Q and ABZ resulted in a lower SASA value than the unbound state (Fig 4E and 4F, dark green and dark pink, respectively). The ligand ABZ and 42Q binding causes the proteins to become more compact and less flexible.

Molecular mechanics generalized born and surface area (MMGBSA) calculations
The MMGBSA technique is widely used to determine the binding energy of ligands to protein molecules. The binding free energy of each of the 3VHE+ABZ and 3VHE+42Q complexes was calculated and the involvement of additional non-bonded interaction energies. The average binding energy of the ABZ bound to 3VHE was found to be -42.07 ± 2.4 kcal/mol, whereas the average binding energy of 42Q was determined to be -80.94 ± 2.25 kcal/mol ( Table 2).

PLOS ONE
G bindCoulomb , G bindCovalent , G bindHbond , G bindLipo , G bindSolvation , and G bindvdW interactions are all examples of non-bonded interactions that affect G bind . G bindvdW , G bindLipo , G bindHbond , and G bindCoulomb energies contributed the most to obtain the average binding energy across all types of interactions. The G bindSolvation and G bindcovalent energies, on the other hand, contributed the least to the final average binding energies. MMGBSA trajectories of 3VHE+ABZ and 3VHE+42Q complexes are further analyzed for positional variances of the respective ligand at the binding pocket. Complex 3VHE with ABZ displayed change in orientation in the ligand position from 0 to 200 ns (Fig 5A). While least arrangement is observed in the 42Q ligand (Fig 5B). MMGBSA trajectory of ABZ and 42Q

PLOS ONE
with 3VHE displayed a very stable configuration throughout the simulation. Moreover, high binding energies were observed due to stable conformation with the protein-ligand complex. Therefore, for the entire simulation studies and MMGBSA analysis, it can be suggested that both ABZ and 42Q have significant affinity for the target 3VHE and can be potential inhibitors of drug molecules.

Analysis of principal component (PCA) with free energy surface/landscape (FEL)
Principal component analysis (PCA) of the MD simulation trajectories for ABZ bound to the VEGFR2 kinase domain was analyzed to interpret the randomized, global motion of the atoms of amino acid residues along with free energy. This analysis interprets the more flexible scattered trajectories due to the protein structure's randomness due to non-correlated global motion. The internal coordinate's mobility into three-dimensional space in the spatial time of 200 ns was recorded in a covariance matrix. The rational motion of each trajectory is interpreted in the form of orthogonal sets or Eigenvectors. MD simulation trajectory of Cα atoms of ABZ bound to VEGFR2 kinase domain displayed more unordered orientation in PC1 and PC2 modes where, the free energies are not converged localized (Fig 6A, blue, as per scale). Following this, PC2 and PC3 mode displayed better order Eigenvalues for the trajectories with partially ordered free energy surface (Fig 6B). The Eigenvectors displayed a less scattered, more correlated and partially converged free energy surface of the last 50 trajectories (Fig 6B). Most of the trajectories finally settled in PC9 and PC10, where the global motion centred toward the origin of the plot (Fig 6C). Eigenvectors are observed. The global movement of MD trajectories toward a positive correlation Eigen vector indicates that the VEGFR2 kinase domain complex is compact and converged. The high periodic global motion was observed along the PC9 and PC10 planes due to the grouping of trajectories in a single free energy surface at the centre of the PCA plot toward origin. Cantering of the trajectories in a single cluster indicates the periodic motion of MD trajectories due to stable conformational global motion (Fig 6C, blue, as per scale).
On the other hand, the 42Q bound VEGFR2 kinase domain displayed a similar observation for PCA as found in ABZ bound state. The MD trajectories in PC1 and PC2 modes displayed more disoriented distribution and less free energy converged zones (Fig 6D, blue), while in PC2 and PC3 modes, trajectories started organizing more into partially definite free energies (Fig 6E, blue). Interestingly increment of PC modes up to PC9 and PC10 displayed more ordered trajectories. The last 50 trajectories (yellow dots) in PC9 and PC10 showed ordered orientation toward the origin of the Eigenvalues (Fig 6F, blue). Ordered orientation signifies the ordered global motion of the trajectories that arose due to a stable converged structure and high free energies toward zero. Therefore, at this stage global minimization achieved.

Conclusion
Albendazole approved by FDA for therapeutic usage as an anthelmintic found to act as VEGFR-2 inhibitor, that represent potential treatment option after investigation. Albendazole binds with conserved residues presents in substrate binding pocket and forms remarkable interactions. Molecular docking and dynamics study are well-established computational techniques that predict drugs' affinity and stability towards the receptor. The current investigation on Albendazole interaction with VEGFR-2 using molecular docking suggested the achievable interaction and possible inhibition. In addition, molecular dynamics simulation suggested greater stability of VEGFR-2 with Albendazole like a standard cocrystal ligand 42Q with plausible and significant binding energies. The results of this study confirm initial reports however additional investigation necessary to check the efficacy of drug by in vitro and in vivo studies followed by clinical trials against ovarian cancer patients. But publishing this information in public domain may help community to fight with ovarian cancer if any potential investigator proceeds for trials. We recommend that this drug candidate be experimentally tested and used as a starting point for further design of a high efficient drug.